Proyecto final ICDS 2020

Consigna:

Realizar una pregunta en su área de interés y responderla con el uso de datos.

Pregunta general:

¿Cuáles son las características de los servicios veterinarios que existen en CABA en relación a la población animal que aspiran satisfacer?

Preguntas secundarias:

¿Cuál es la distribución de los servicios?
¿Qué tipos de servicios veterinarios brindan los establecimientos en CABA?
¿La financiación es mayoritariamente pública o privada?
¿Qué relación hay entre la distribución de la población animal y los servicios?
¿Hay zonas de CABA donde coexistan servicios públicos de diferentes reparticiones?
¿Los servicios públicos están donde hacen falta?

MARCO TEORICO:

Diversos factores han favorecido el incremento del número de animales de compañía en las grandes ciudades. El veterinario es el profesional llamado a participar activamente en el complejo entramado de la relación entre propietario y animal. Su participación debe involucrar los aspectos de sanidad animal, la nutrición y la educación a propietarios sobre la tenencia y manejo del animal compañía.
Según proyecciones de la Dirección General de Estadística y Censos (DGEyC) para el año 2019 había un estimado de 3.072.029 de habitantes en el área de 203 kilómetros cuadrados de la CABA. De acuerdo con los resultados de la Encuesta Anual de Hogares (EAH) 2018, se estima que existen cerca de 475.000 perros y alrededor de 295.000 gatos en los hogares de la CABA, ésta encuesta excluyó del análisis a los animales callejeros sin dueño. La relación entre la cantidad de animales y la cantidad de habitantes presenta diferencias entre las comunas.
En la ciudad existen 2 servicios veterinarios públicos gratuitos: el Instituto de Zoonosis Luis Pasteur (IZLP) dependiente del Ministerio de Salud del Gobierno de la Ciudad Autónoma de Buenos Aires y Mascotas de la Ciudad de la Agencia de Protección Ambiental dependiente de la Secretaría de Ambiente. Como medidas sobre la población animal se enfocan en la vacunación antirrábica obligatoria por ley, desparasitación, esterilizaciones quirúrgicas para control poblacional y atención integral, con foco en detección de zoonosis en el caso del IZLP.
El hospital escuela de la Facultad de Ciencias Veterinaria de la Universidad de Buenos Aires (FCV-UBA) presenta diferencias: se trata de un servicio arancelado a menor costo que en los locales privados, incluye todas las especialidades médicas y su principal objetivo es ser una herramienta de enseñanza que presta un servicio a la comunidad. Para atender exigen que los animales presenten la vacuna antirrábica y no realizan ovariectomias/orquiectomias para control poblacional, sino sólo sí el cuadro clínico del paciente lo amerita.
Con respecto a los privados, el estudio más exhaustivo hasta la fecha fue realizado por el cuerpo de inspectores del Consejo Profesional de Médicos Veterinarios (CPMV) que durante 2019 contabilizó 1133 establecimientos en actividad. De los cuales aproximadamente la mitad brinda atención clínica y unos 14 efectores brindan atención las 24hs. La OPS usa un indicador de disponibilidad para medicina humana llamado “Cantidad de establecimientos de salud por cada mil habitantes” que se calcula como “Número de establecimientos de salud en una determinada área geográfica y/o dependencia administrativa para un periodo dado, en un determinado país, por cada 1000 habitantes.”
Extrapolando este indicador a nuestro campo podríamos calcular algo similar como una “Razón de Establecimientos Veterinarios”:

1. En el numerador corresponde un conteo del número de establecimientos, solo se tendrán en cuenta aquellos que presten atención clínica por parte de un veterinario, en el caso de los publicos se contaran aquellos puntos con una sede fija con intraestrucura (quedan excluidos los puestos móviles) y en cuanto a los privados se definiran siuiendo la lista de categorías del CPMV
2. En el denominador se calculará sobre el estimado de animales por comuna. Para obtener dicha estimación cruzaremos los datos de población por comuna año 2019 con los obtenidos por la EAH 2018, dividido por 1000.

FÓRMULA:
“Número de establecimientos veterinarios en comuna”x" de la CABA periodo 2019 / (Número estimado de animales en comuna “x” / 1000)”

Dada la variabilidad de características y factores a considerar, el uso de herramientas como los Sistemas de Información Geográfica (SIG) son altamente eficaces por su versatilidad para incorporar información siendo una opción adecuada para estudiar la distribución en relación al espacio de los distintos prestadores de acuerdo a su ubicación en las poblaciones que aspiran satisfacer.

Mascotas de la Ciudad
Instituto de Zoonosis Luis Pasteur
Facultad de Ciencias Veterinarias UBA
Consejo Profesional de Medicos Veterinarios

---

Librerias usadas:

library(tidyverse)
library(rmarkdown)
library(RUMBA)
library(sf)
library(readxl)
library(ggmap)
library(osmdata)
library(writexl)
library(leaflet)

Bases:

Comenzamos cargando las bases, 3 formato excell y 1 geojson conteniendo las comunas de CABA:

#servet <- read_excel("C:/Users/Federico/Downloads/servet.xlsx") esta es la base inicial

servet <- read_excel("C:/Users/Federico/Downloads/servet_CONCOORDENADAS.xlsx") #esta es la base final

comunas <- st_read("C:/Users/Federico/Downloads/CABA_comunas.geojson") #base con gatos geo

## Reading layer `comunas' from data source `C:\Users\Federico\Downloads\CABA_comunas.geojson' using driver `GeoJSON'
## Simple feature collection with 15 features and 6 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -58.53152 ymin: -34.70529 xmax: -58.33515 ymax: -34.52649
## geographic CRS: WGS 84

datos_comunas <- read_excel("C:/Users/Federico/Downloads/poblacion_animal_comuna.xlsx") #poblacion animal

datos_eah <- read_excel("C:/Users/Federico/Downloads/resumen_eah_mascotas.xlsx") #indices de interes de la EAH

#Exploramos
head(servet)

## # A tibble: 6 x 13
##      id unidad direc barrio comuna zona  origen_financia~ movil servicio
##   <dbl> <chr>  <chr> <chr>   <dbl> <chr> <chr>            <dbl> <chr>   
## 1     1 Masco~ Av. ~ Puert~      1 CENT~ publico              0 At. Cli~
## 2     2 Masco~ Av. ~ <NA>        1 CENT~ publico              1 At. Cli~
## 3     3 Masco~ Chil~ <NA>        1 CENT~ publico              1 At. Cli~
## 4     4 Masco~ Av. ~ <NA>        1 CENT~ publico              1 At. Cli~
## 5     5 Masco~ Para~ <NA>        2 NORTE publico              1 At. Cli~
## 6     6 Masco~ Dr. ~ <NA>        2 NORTE publico              1 At. Cli~
## # ... with 4 more variables: periodicidad <chr>, address_normalised <chr>,
## #   lon <dbl>, lat <dbl>

head(comunas)

## Simple feature collection with 6 features and 6 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -58.53152 ymin: -34.67451 xmax: -58.3707 ymax: -34.56827
## geographic CRS: WGS 84
##                                                                           BARRIOS
## 1                                                                        RECOLETA
## 2                                                                 ALMAGRO - BOEDO
## 3                                                                       CABALLITO
## 4                                                       FLORES - PARQUE CHACABUCO
## 5                                         LINIERS - MATADEROS - PARQUE AVELLANEDA
## 6 FLORESTA - MONTE CASTRO - VELEZ SARSFIELD - VERSALLES - VILLA LURO - VILLA REAL
##   PERIMETRO     AREA COMUNAS ID         OBJETO                       geometry
## 1  21452.84  6317265       2  1 LIMITE COMUNAL MULTIPOLYGON (((-58.38 -34....
## 2  12323.43  6660603       5  2 LIMITE COMUNAL MULTIPOLYGON (((-58.41287 -...
## 3  10990.96  6851029       6  3 LIMITE COMUNAL MULTIPOLYGON (((-58.43061 -...
## 4  17972.26 12422901       7  4 LIMITE COMUNAL MULTIPOLYGON (((-58.452 -34...
## 5  21411.74 16505306       9  5 LIMITE COMUNAL MULTIPOLYGON (((-58.51925 -...
## 6  18332.04 12656557      10  6 LIMITE COMUNAL MULTIPOLYGON (((-58.48834 -...

head(datos_comunas)

## # A tibble: 6 x 5
##   comuna poblacion pc100hab gc100hab porc_NBI
##    <dbl>     <dbl>    <dbl>    <dbl> <chr>   
## 1      1    255457     11.5      8   15.9    
## 2      2    149510      9        6.5 2       
## 3      3    193115     12        9.5 11.9    
## 4      4    239712     21.5      9.7 12.7    
## 5      5    187348     14.5     11   6.1     
## 6      6    185271     11.5      7.5 2.2

head(datos_eah)

## # A tibble: 6 x 27
##   comuna perro_nodespara~ gato_nodesparas~ total_gato_desp~ total_perro_des~
##    <dbl>            <dbl>            <dbl>            <dbl>            <dbl>
## 1      1              132              148              247              326
## 2      2               40               68              105              121
## 3      3              124               92              178              267
## 4      4              260              196              344              593
## 5      5              100               88              163              211
## 6      6               84               64              117              200
## # ... with 22 more variables: perro_vacunado <dbl>, perro_novacunado <dbl>,
## #   total_perro_vacunacion <dbl>, gato_vacunado <dbl>, gato_novacunado <dbl>,
## #   total_gato_vacunacion <dbl>, perro_castrado <dbl>, perro_nocastrado <dbl>,
## #   gato_castrado <dbl>, gato_nocastrado <dbl>, total_perros_castracion <dbl>,
## #   total_gatos_castracion <dbl>, cuida_animalcomunitario <dbl>,
## #   nocuida_animalcomunitario <dbl>, total_cuidado_animalcomunitario <dbl>,
## #   porc_cuida_animalcomunitario <dbl>, adq_compro <dbl>, adq_regalo <dbl>,
## #   adq_recogido <dbl>, adq_adoptado <dbl>, adq_criapropia <dbl>,
## #   adq_total <dbl>

Tenemos una base de 983 datos correspondientes a lugares en CABA donde se prestan servicios veterinarios.
Ademas otras bases donde se resume la informacion obtenida con respecto a las mascotas de la ciudad, podemos estimar la poblacion animal, el porcentaje de animales castrados, los que no fueron vacunados con la antirrabica ni desparasitados en los ulitmos 12 meses al momento de la encuesta para las 15 comunas de la CABA.

Transformacion:

Antes de ponernos artísticos con los mapas, necesitamos hacer unas transformaciones de los datos tal y como nos los dieron para que sean más fácil de entender y manipular. Veamos como estaría distribuida la población animal en la CABA. Lo primero es hacer un join entre los shapes de las comunas y los datos que les corresponden, extraídos de la Encuesta Anual de Hogares 2018 - Modulo de Tenencia Responsable de animales de compañía:

comunas <- left_join(comunas,datos_comunas,by = c("COMUNAS" = "comuna"))

#Ahora calculamos en base a la poblacion humana, la poblacion animal. Respeto la separacion entre caninos y felinos porque hay otra base de interes que contiene datos diferenciados entre estas dos especies de animales de compañia:

comunas <- comunas %>% 
  mutate(caninos_comunas = round((poblacion/100) * pc100hab),
         felinos_comunas = round((poblacion/100) * gc100hab),
         animales_comunas = round(caninos_comunas + felinos_comunas)) #round = para manejar enteros, ya que hablamos de individuos!

#Para chequear que los calculos sean correctos, hacemos unas sumatorias.
sum(comunas$caninos_comunas)

## [1] 475215

sum(comunas$felinos_comunas)

## [1] 297429

sum(comunas$animales_comunas)

## [1] 772644

#Según la E.A.H. hay alrededor de 475000 caninos y 300000 felinos en los hogares de CABA. Vamos bien!
head(comunas)

## Simple feature collection with 6 features and 13 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -58.53152 ymin: -34.67451 xmax: -58.3707 ymax: -34.56827
## geographic CRS: WGS 84
##                                                                           BARRIOS
## 1                                                                        RECOLETA
## 2                                                                 ALMAGRO - BOEDO
## 3                                                                       CABALLITO
## 4                                                       FLORES - PARQUE CHACABUCO
## 5                                         LINIERS - MATADEROS - PARQUE AVELLANEDA
## 6 FLORESTA - MONTE CASTRO - VELEZ SARSFIELD - VERSALLES - VILLA LURO - VILLA REAL
##   PERIMETRO     AREA COMUNAS ID         OBJETO poblacion pc100hab gc100hab
## 1  21452.84  6317265       2  1 LIMITE COMUNAL    149510      9.0      6.5
## 2  12323.43  6660603       5  2 LIMITE COMUNAL    187348     14.5     11.0
## 3  10990.96  6851029       6  3 LIMITE COMUNAL    185271     11.5      7.5
## 4  17972.26 12422901       7  4 LIMITE COMUNAL    241484     16.0      9.0
## 5  21411.74 16505306       9  5 LIMITE COMUNAL    171062     22.0     12.0
## 6  18332.04 12656557      10  6 LIMITE COMUNAL    170497     18.0      8.5
##   porc_NBI                       geometry caninos_comunas felinos_comunas
## 1        2 MULTIPOLYGON (((-58.38 -34....           13456            9718
## 2      6.1 MULTIPOLYGON (((-58.41287 -...           27165           20608
## 3      2.2 MULTIPOLYGON (((-58.43061 -...           21306           13895
## 4      8.6 MULTIPOLYGON (((-58.452 -34...           38637           21734
## 5      4.2 MULTIPOLYGON (((-58.51925 -...           37634           20527
## 6      3.5 MULTIPOLYGON (((-58.48834 -...           30689           14492
##   animales_comunas
## 1            23174
## 2            47773
## 3            35201
## 4            60371
## 5            58161
## 6            45181

Segundo, vamos a completar la base de los servicios con las zonas en las que la E.A.H. divide a la CABA:

CENTRO <- c("1", "3", "5", "6", "7", "11", "12", "15")

NORTE <- c("2", "13", "14")

SUR <- c("4", "8", "9", "10")

#servet <- servet %>%
#  mutate(zona = case_when(comuna %in% CENTRO ~ "CENTRO",
#                          comuna %in% NORTE ~ "NORTE",
#                          comuna %in% SUR ~ "SUR"))
head(servet)

## # A tibble: 6 x 13
##      id unidad direc barrio comuna zona  origen_financia~ movil servicio
##   <dbl> <chr>  <chr> <chr>   <dbl> <chr> <chr>            <dbl> <chr>   
## 1     1 Masco~ Av. ~ Puert~      1 CENT~ publico              0 At. Cli~
## 2     2 Masco~ Av. ~ <NA>        1 CENT~ publico              1 At. Cli~
## 3     3 Masco~ Chil~ <NA>        1 CENT~ publico              1 At. Cli~
## 4     4 Masco~ Av. ~ <NA>        1 CENT~ publico              1 At. Cli~
## 5     5 Masco~ Para~ <NA>        2 NORTE publico              1 At. Cli~
## 6     6 Masco~ Dr. ~ <NA>        2 NORTE publico              1 At. Cli~
## # ... with 4 more variables: periodicidad <chr>, address_normalised <chr>,
## #   lon <dbl>, lat <dbl>

Esta diferenciación se adapta a los relevamientos de la EAH ya que tienen distintos porcentajes en las dimensiones del cuidado de las mascotas, proporción de hogares que colaboran con el cuidado de animales callejeros, ratios de perros por persona, atención veterinaria de los animales de compañía y distintos niveles de ingresos que configuran escenarios heterogéneos que deben ser ponderados en el análisis final.

Esta división en principio es netamente territorial, pero por los índices se aprecian diferencias entre los comportamientos de los propietarios entre las zonas. Mande mails para tratar de profundizar esta información y ver si puedo rescatar los datos de forma más desagregada, idealmente llegar hasta Barrio. Me contestaron que la metodología lo impide, consulte por datos a nivel comuna y estoy a la espera de los mismos.

Transformamos la columna dirección en datos longitud y latitud:

#servet <- mutate_USIG_geocode(servet, "direc")

#head(servet)
#summary(servet)

En este trabajo final ya están depurados, pero la primera inspección de los sitios georreferenciados, revelaronn que existe una proporción significativa de registros para los cuales no fue posible traducir la dirección disponible en coordenadas precisas. En algunos casos estaban mal las direcciones (cargadas a mano la mayoria) y en otros aun chequeandolo por google no me las reconocian.

En cuanto a los servicios púbicos también se encuentran algunas posiciones duplicadas (Puntos de atención georreferenciados con las mismas coordinadas pero distinto servicio) que corresponden a lugares frecuentes donde un mismo servicio va a realizar actividades distintas, o donde dos reparticiones prestan atención en momentos distintos, esto es importante porque forma parte de los objetivos de análisis.

Tuve que realizar una limpieza manual de las direcciones para mejorar su precisión y volver a geolocalizarlas. Fue bastante laborioso, pero me asegure de que se aprovecharon el máximo de datos que me fueron proporcionados. Este paquete de SIG para CABA demostró tener varios baches para geolocalizar las direcciones cuando caían en Barrios Vulnerables, particularmente tuve numerosos “NA” en Soldati y Retiro. Para no perder los datos buscaba y reemplazaba la dirección la esquina más cercana que me reconociera el paquete, elegí perder precisión, pero conservar la mayor cantidad de datos, asegurándome de que no caigan en otro barrio/comuna y me distorsionen el análisis.

Aplicamos la función para transformar los datos en objetos geométricos, pero creando otro dataset:

servet_geo <- servet %>% 
  filter(!is.na(lon), !is.na(lat)) %>% #filtrado de valores ausentes
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

head(servet) #se crearon 3 nuevas columnas, una con la direccion normalizada y las restantes con los datos de "lon" y "lat"

## # A tibble: 6 x 13
##      id unidad direc barrio comuna zona  origen_financia~ movil servicio
##   <dbl> <chr>  <chr> <chr>   <dbl> <chr> <chr>            <dbl> <chr>   
## 1     1 Masco~ Av. ~ Puert~      1 CENT~ publico              0 At. Cli~
## 2     2 Masco~ Av. ~ <NA>        1 CENT~ publico              1 At. Cli~
## 3     3 Masco~ Chil~ <NA>        1 CENT~ publico              1 At. Cli~
## 4     4 Masco~ Av. ~ <NA>        1 CENT~ publico              1 At. Cli~
## 5     5 Masco~ Para~ <NA>        2 NORTE publico              1 At. Cli~
## 6     6 Masco~ Dr. ~ <NA>        2 NORTE publico              1 At. Cli~
## # ... with 4 more variables: periodicidad <chr>, address_normalised <chr>,
## #   lon <dbl>, lat <dbl>

class(servet)

## [1] "tbl_df"     "tbl"        "data.frame"

head(servet_geo) #transforma los datos "lon" y "lat" en "geometry"

## Simple feature collection with 6 features and 11 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -58.39493 ymin: -34.62212 xmax: -58.35644 ymax: -34.58339
## geographic CRS: WGS 84
## # A tibble: 6 x 12
##      id unidad direc barrio comuna zona  origen_financia~ movil servicio
##   <dbl> <chr>  <chr> <chr>   <dbl> <chr> <chr>            <dbl> <chr>   
## 1     1 Masco~ Av. ~ Puert~      1 CENT~ publico              0 At. Cli~
## 2     2 Masco~ Av. ~ <NA>        1 CENT~ publico              1 At. Cli~
## 3     3 Masco~ Chil~ <NA>        1 CENT~ publico              1 At. Cli~
## 4     4 Masco~ Av. ~ <NA>        1 CENT~ publico              1 At. Cli~
## 5     5 Masco~ Para~ <NA>        2 NORTE publico              1 At. Cli~
## 6     6 Masco~ Dr. ~ <NA>        2 NORTE publico              1 At. Cli~
## # ... with 3 more variables: periodicidad <chr>, address_normalised <chr>,
## #   geometry <POINT [°]>

class(servet_geo)

## [1] "sf"         "tbl_df"     "tbl"        "data.frame"

Me interesa tener 2 datasets con la misma información, pero uno clase “data.frame” y otro clase “sf, data.frame” porque al graficar voy a usar diferentes paquetes que toman uno u otro objeto, para no andar convirtiendo y reconvirtiendo los dejo duplicados pero diferenciados con el sufijo "_geo".

Filtros:

Acá me propuse filtrarlos por diversas características que me parecían relevantes para luego poder graficarlas o cruzarlas con otros datos:

#Separando las privadas
privadas <- servet %>% 
  filter(origen_financiamiento == "privada")

privadas_geo <- privadas %>%
  filter(!is.na(lon), !is.na(lat)) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

#Dentro de las privadas hay un subconjunto interesante: las 24hs. Son pocas comparadas al resto, unas 14.

privadas_24hs <- privadas %>% 
  filter(periodicidad == "24hs")
head(privadas_24hs)

## # A tibble: 6 x 13
##      id unidad direc barrio comuna zona  origen_financia~ movil servicio
##   <dbl> <chr>  <chr> <chr>   <dbl> <chr> <chr>            <dbl> <chr>   
## 1   199 Efect~ FREN~ <NA>        2 NORTE privada              0 clinica~
## 2   332 Efect~ BOED~ <NA>        5 CENT~ privada              0 clinica~
## 3   399 Efect~ AV. ~ <NA>        6 CENT~ privada              0 clinica~
## 4   401 Efect~ DIRE~ <NA>        6 CENT~ privada              0 clinica~
## 5   469 Efect~ EVA ~ <NA>        7 CENT~ privada              0 clinica~
## 6   509 Efect~ MART~ <NA>        8 SUR   privada              0 clinica~
## # ... with 4 more variables: periodicidad <chr>, address_normalised <chr>,
## #   lon <dbl>, lat <dbl>

privadas_24hs_geo <- privadas_24hs %>%
  filter(!is.na(lon), !is.na(lat)) %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

#Veamos que servicios declaran los prestadores
unique(servet$servicio)

##  [1] "At. Clinica + Castracion"                                                                           
##  [2] "At. Clinica + Castracion + Vac. Antirrabica"                                                        
##  [3] "Vac. Antirrabica"                                                                                   
##  [4] "Desparasitacion + Vac. Antirrabica"                                                                 
##  [5] "Castracion"                                                                                         
##  [6] "At. Clinica + Castracion+ Desparasitacion + Diag. Laboratorio + Toma de Muestras + Vac. Antirrabica"
##  [7] "At. Clinica + Desparasitacion + Toma de Muestras + Diag. Laboratorio"                               
##  [8] "At. Clinica + Desparasitacion + Toma de Muestras + Vac. Antirrabica"                                
##  [9] "centro + zooterapicos"                                                                              
## [10] "consultorio + petshop + zooterapicos"                                                               
## [11] "otros"                                                                                              
## [12] "petshop"                                                                                            
## [13] "petshop + zooterapicos"                                                                             
## [14] "acuario"                                                                                            
## [15] "clinica_internacion + zooterapicos"                                                                 
## [16] "clinica_simple + zooterapicos"                                                                      
## [17] "consultorio"                                                                                        
## [18] "laboratorio"                                                                                        
## [19] "centro + petshop + zooterapicos"

#Para confeccionar la Razón de Establecimiento Veterinarios tuve que definir 2 parámetros: pirmero que sean establecimientos fijos, que cuenten con infraestructuras, no puestos móviles; segundo que presten atención clínica. En base a estas dos características filtro mi base:

establecimientos <- servet %>% 
  filter(movil == "0", servicio  %in% c("At. Clinica + Castracion", 
                                        "At. Clinica + Castracion + Vac. Antirrabica",
                                        "At. Clinica + Castracion+ Desparasitacion + Diag. Laboratorio + Toma de Muestras + Vac. Antirrabica", 
                                        "At. Clinica + Desparasitacion + Toma de Muestras + Diag. Laboratorio",
                                        "centro + zooterapicos",
                                        "consultorio + petshop + zooterapicos",
                                        "clinica_internacion + zooterapicos",
                                        "clinica_simple + zooterapicos",
                                        "consultorio",
                                        "centro + petshop + zooterapicos"))

establecimientos_geo <- establecimientos %>%
  filter(!is.na(lon), !is.na(lat)) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

#Vamos con los publicos
publicos <- servet %>% 
  filter(origen_financiamiento == "publico")

#Me interesa de alguna forma resaltar aquellos puntos que tengan más frecuencia de prestación en el servicio. Veamos que declaran: 
unique(servet$periodicidad)

##  [1] "lunes a viernes y fin de semana por medio"
##  [2] "mensual"                                  
##  [3] "bimensual"                                
##  [4] "trimestral"                               
##  [5] "quincenal"                                
##  [6] "cuatrimestral"                            
##  [7] "semestral"                                
##  [8] "anual"                                    
##  [9] "bianual"                                  
## [10] "semanal"                                  
## [11] "dias hábiles"                             
## [12] "horario comercial variable"               
## [13] "24hs"

#Arbitrariamente definí como presencia significativa que sea atención mínimamente mensual en dicho punto: 

puestos_frecuentes <- servet %>% 
  filter( periodicidad %in% c("lunes a viernes y fin de semana por medio", 
                              "dias hábiles", 
                              "semanal", 
                              "quincenal", 
                              "mensual"))

puestos_frecuentes_geo <- puestos_frecuentes %>%
  filter(!is.na(lon), !is.na(lat)) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

Primer vistazo a los datos: MAPA BASICO DE CABA:

A esta altura tenemos una cantidad importante de elementos, asi que empecemos con unos mapas simples para tener nociones:

#Obtenemos la capa base de la Ciudad Autonoma de Buenos Aires:

bbox <- getbb("Ciudad Autónoma de Buenos Aires, Argentina") #especificamos ciudad

CABA_blanco <- get_stamenmap(bbox = bbox,
                      maptype = "toner-hybrid",zoom=12) #especificamos tipo de mapa, este al ser un blanco me queda mas comodo para ubicar las nubes de densidad

CABA_terrain <- get_stamenmap(bbox = bbox,
                      maptype = "terrain",zoom=12) #este me gusta mas para la visualizacion de puntos ya que tiene un adecuado balance de vias importantes y zonas verdes.

ggmap(CABA_blanco)

ggmap(CABA_terrain)

#Empecemos viendo la totalidad de los puntos para tener una primera impresión de los datos:
ggmap(CABA_terrain) +
  geom_sf(data=servet_geo, 
          color="blue", 
          size=1, 
          inherit.aes=FALSE)+
  labs( x = "Longitud", 
        y = "Latitud",
        title = "Servicios veterinarios en CABA",
        caption = "Fuente: CPMV, IZLP, Mascotas Ciudad, GoogleMaps")

Se ve que están distribuidos por toda la ciudad, me arriesgo a decir que se nota una menor densidad en la zona sur…

Mapa de poblacion animal en CABA:

Antes de sumergirnos en los mapas de los servicios, veamos como se distribuye la población animal en la CABA:

#Primero una tabla que da pie al mapa para tener los numeros presentes:

ggplot(comunas) + 
  geom_bar(aes(x= reorder(COMUNAS, animales_comunas), 
               weight=animales_comunas)) + 
  labs( x = "Comuna", 
        y = "Población animal estimada 2019",
        title = "Cantidad estimada de animales por comuna",
        caption = "Fuente: Informe módulo de Tenencia responsable y sanidad de perros y gatos 2018") +
  coord_flip()+
  theme_minimal()

#El mapa
ggplot() + 
  geom_sf(data=comunas, aes(fill= animales_comunas), color = "black") +
  scale_fill_viridis_c() +
    theme_minimal()+
   labs(x = "Longitud", 
        y = "Latitud",
       title ="Cantidad estimada de animales por comuna",  
       subtitle = "Ciudad Autónoma de Buenos Aires", 
       fill = "Numero estimado de animales",
       caption = "Fuente: Informe módulo de Tenencia responsable y sanidad de perros y gatos 2018")

En la comuna 4(sur) tenemos un estimado de 75.000 animales, seguido por comuna 8(sur) con 65.000, le siguen comuna 7(centro), 9(sur), 12(centro) y 13(norte) con alrededor de 60.000 cada una.

Facetando:

Veamos por prestadores publicos (Mascotas de la Ciudad, Instituto Pasteur y la UBA) versus los prestadores privados sin distincion de rubro:

ggmap(CABA_terrain) +
    geom_point(data = servet, 
               aes(x = lon, y = lat), 
               color="blue", 
               size=1, inherit.aes=FALSE) + 
  facet_wrap(~ unidad)+
  labs(x = "Longitud", 
       y = "Latitud",
       title ="Servicios veterinarios 2019",  
       subtitle = "Ciudad Autónoma de Buenos Aires",
       caption = "Fuente: CPMV, IZLP, Mascotas Ciudad, GoogleMaps")

#Aca claramente comienza una divergencia entre los prestadores. Ademas de la distribucion, algo evidente es la cantidad.

Estos detalles se puedeen apreciar mejor con otro tipo de grafico:

Grafico heatmap o nube de densidad:

ggmap(CABA_blanco) + 
  stat_density2d(data = servet, 
                 aes(x= lon, y= lat, fill=stat(level)), 
                 alpha=0.4, 
                 geom="polygon")+
  scale_fill_viridis_c()+
  facet_wrap(~ origen_financiamiento)+
  labs(x = "Longitud", 
       y = "Latitud",
       title = "Densidad de servicios veterinarios CABA",
       subtitle = "Dividido por origen de financiamiento",
       caption = "Fuente: CPMV, IZLP, Mascotas Ciudad, GoogleMaps")

Acá vemos una clarísima diferencia en la distribución de ambos servicios.

Los privados están ampliamente distribuidos en las zonas norte y centro con menor presencia en la zona sur. Hay 2 puntos de alta concentración en las zonas de las comunas 2(norte), 3(centro), 5(centro), 13(norte) y 14(norte).

Los públicos cubren casi toda la CABA con 2 puntos de alta concentración en las comunas 4(sur) y 8(sur).

ggmap(CABA_blanco) + #mismos datos que el ultimo grafico, pero como una nube de densidad de puntos,
  stat_density2d(data = servet, #este dataset tiene en cuenta TODOS los efectores, fijos y moviles
                 aes(x= lon, y= lat, fill=stat(level)), #y TODOS los servicios ofrecidos
                 alpha=0.4, 
                 geom="polygon")+
  scale_fill_viridis_c()+
  facet_wrap(~ unidad)+
  labs(x = "Longitud", 
       y = "Latitud",
       title = "Densidad de servicios veterinarios CABA",
       subtitle = "Dividido por prestadores")

#El hospital escuela sale vacio por ser un unico punto, lo podemos evidenciar mostrando los puntos:

ggmap(CABA_blanco) +
  stat_density2d(geom = "polygon", 
                 data = servet, 
                 aes(x= lon, 
                     y= lat, 
                     fill=stat(level)),
                 alpha=0.4)+
  geom_point(data = servet, 
               aes(x = lon, 
                   y = lat), 
               color="red", 
               size=1,
             alpha=0.25,
             inherit.aes=FALSE)+
  scale_fill_viridis_c()+
  facet_wrap(~ unidad)+
  labs(x = "Longitud", 
       y = "Latitud",
       title = "Densidad de servicios veterinarios CABA",
       subtitle = "Dividido por prestadores, cada punto rojo es un lugar de atención")

#Último donde pongo de manifiesta los puestos FIJOS vs MOVILES en los públicos:
ggmap(CABA_blanco) +
  stat_density2d(geom = "polygon", #lo mismo
                 data = servet, 
                 aes(x= lon, 
                     y= lat, 
                     fill=stat(level)),
                 alpha=0.25)+
  geom_point(data = establecimientos, #pero con puntos que reflejan PUESTOS FIJOS que prestan ATENCION CLINICA
               aes(x = lon, 
                   y = lat), 
               color="red",
               size=2,
             alpha=0.25,
             inherit.aes=FALSE)+
  scale_fill_viridis_c() +
  facet_wrap(~ unidad)+
  labs(x = "Longitud", 
       y = "Latitud",
       title = "Densidad de servicios veterinarios CABA",
       subtitle = "Dividido por prestadores, puntos rojos son lugares fijos de atencion")

Mascotas de la Ciudad afirma cubrir las 15 comunas de CABA y visualmente vemos un patrón homogéneo, están presentes en toda la ciudad salvo por dos áreas blancas al norte. Tengamos en cuenta que estos datos corresponden a un año, las posiciones suelen rotar yyyyy un factor importantísimo (a mi criterio) que el 2019 fue un año electoral.

El IZLP en sus orígenes estuvo concentrado en su sede central de Parque Centenario, luego la División de Acciones Comunitarias de la Salud (DACS o el Área Programática de los hospitales) empezó a ir a ciertos puntos de la ciudad con puestos móviles, finalmente se sumaron las acciones de la Residencia de Veterinaria en Salud Publica que atienden semanalmente en los Cesacs 24 y 43. Tiene un patrón más heterogéneo, concentrado en el sur y hay varios puntos desperdigados por el resto de CABA.

La Facultad de Ciencias Veterinarias de la UBA en el barrio de agronomía fue detallada en el marco teórico.

Mientras que las privadas vemos una altísima concentración en la zona norte y centro, un poco menos de densidad hacia el oeste y una zona más vacía en el sur. Esto puede ser medio engañoso porque hay muchas zonas de parques (véase la mancha verde del cementerio de la chacarita), veremos que dicen cuando los sometemos a otros gráficos.

Si bien vimos que los públicos tienen buena distribución sobre el territorio, veamos qué pasa si mostramos los puntos de ATENCION FRECUENTE que definimos y filtramos, sobre la población animal estimada por comunas:

ggplot() + 
  geom_sf(data=comunas, aes(fill= animales_comunas), color = "black") +
  scale_fill_viridis_c() +
    theme_minimal()+
  geom_sf(data = puestos_frecuentes_geo, color = "red", size = 3)+ # puntos de atencion frecuente
  labs(x = "Longitud", 
       y = "Latitud",
       title ="Cantidad estimada de animales por comuna",  
       subtitle = "Ciudad Autónoma de Buenos Aires", 
       fill = "Numero estimado de animales",
       caption = "Fuente: Informe módulo de Tenencia responsable y sanidad de perros y gatos 2018")

Estan distribuidos por gran parte de CABA, pero no sabemos si son suficientes para la demanda de la población.

Establecimientos privados:

#Según los servicios que declaran: 

ggmap(CABA_blanco) +
  geom_sf(data=privadas_geo, 
          color="blue", 
          size=1, 
          inherit.aes=FALSE) +
  facet_wrap(~ servicio)+
  labs(x = "Longitud", 
       y = "Latitud",
       title ="Servicios según prestaciones que declara",  
       subtitle = "Ciudad Autónoma de Buenos Aires", 
       caption = "Fuente: Consejo Profesional de Médicos Veterinarios CABA")

Claramente hay una abundancia de consultorios con petshops mas venta de zooterapicos que supera a todas las otras categorias, seguido de petshop + zooterapicos.

#Veamos solo los que prestan atencion clinica:
ggmap(CABA_blanco) +
  geom_sf(data=establecimientos_geo, 
          color="blue", 
          size=1, 
          inherit.aes=FALSE)+
  labs(x = "Longitud", 
       y = "Latitud",
       title ="Locales que prestan atención veterinaria clínica",  
       subtitle = "Ciudad Autónoma de Buenos Aires",
       caption = "Fuente: CPMV CABA")

Están ampliamente distribuidos, con mayor densidad centro/norte y tienden a estar sobre vías principales.

Calculo de Razón de Establecimientos Veterinarios:

Es el índice que adoptamos para medir disponibilidad.

establecimientos_comuna <- establecimientos %>% #solo se contemplan los establecimientos FIJOS que hacen atención clínica
    group_by(comuna) %>% #agrupamos por comuna
    summarize(establecimientos_con_clinica = n()) #sumatoria de los elementos del grupo

comunas <- left_join(comunas,establecimientos_comuna,by = c("COMUNAS" = "comuna")) # lo incorporamos a nuestro dataframe comuna

comunas <- comunas %>% 
  mutate(REV = establecimientos_con_clinica/ (animales_comunas/1000)) #fórmula

ggplot() + 
  geom_sf(data=comunas, aes(fill= REV), color = "black") +
  scale_fill_viridis_c() +
    theme_minimal()+
  labs(x = "Longitud", 
       y = "Latitud",
       title ="Razón de Establecimientos Veterinarios",  
       subtitle = "REV = Establecimientos con atención clinica / (Población animal comuna/1000",
       caption = "Fuente: Elaboración personal. Federico Krapp 2020")

Los servicios en relación a la población están concentrados en las comunas 2(norte), 6(centro) y 14(norte).

Las 2 comunas que más población animal tienen, son las 2 que tienen bajísima oferta de servicios veterinarios, que recordemos son alrededor de 600 privados y 4 publicos.

¿Qué factor habrá configurado esto? ¿Las zonas con mayor oferta de servicios, que suelen incluir esterilizaciones, han ejercido presión sobre las poblaciones controlando su número? Si vemos la EAH el porcetanje de animales castrados en CABA ronda el 50%, asi que descartado.

¿Qué otros factores predicen la localización de los servicios veterinarios que puedan explicar esta concentración?

Uno importante que se me ocurre es el poder adquisitivo de la zona, una característica dentro de nuestra profesión es que más allá de la vocación de servicio hacia los animales y la calidad profesional, una buena administración comercial es crucial para el desarrollo (y en algunos casos, la supervivencia) de la oferta de servicios. En ese sentido, una realidad medio chocante es que somos un bien de lujo, porque a diferencia de la salud humana que está garantizada por el Estado y tiene fuerte presencia en CABA, los servicios veterinarios públicos no tienen tanto desarrollo y relevancia, siendo que otros países están mucho mejor articulados con la salud humana bajo el paradigma de Una Salud. Sigamos.

Puestos publicos

publicos_comuna <- publicos %>% #en este caso quise tambien ver donde se concentraron los publicos por comuna en el 2019
  group_by(comuna) %>% 
  summarize(puestos_publicos_anio = n())

publicos_comuna

## # A tibble: 15 x 2
##    comuna puestos_publicos_anio
##     <dbl>                 <int>
##  1      1                     5
##  2      2                     4
##  3      3                     1
##  4      4                    32
##  5      5                     2
##  6      6                     6
##  7      7                    19
##  8      8                    44
##  9      9                    11
## 10     10                     8
## 11     11                     3
## 12     12                     9
## 13     13                     4
## 14     14                     5
## 15     15                     7

comunas <- left_join(comunas,publicos_comuna,by = c("COMUNAS" = "comuna"))

#Veamos un resumen
ggplot(comunas) + 
  geom_bar(aes(x= reorder(COMUNAS, puestos_publicos_anio), 
               weight=puestos_publicos_anio)) +
  theme_minimal()+
  labs(x = "Comunas",
       y ="Puestos móviles realizados en el 2019",
       title ="Cantidad de intervenciones con puestos móviles",
       caption = "Fuente: IZLP, Mascotas de la ciudad")

#Ahora el mapa
ggplot() + 
  geom_sf(data=comunas, aes(fill= puestos_publicos_anio), color = "black") +
  scale_fill_viridis_c() +
    theme_minimal()+
   labs(x = "Longitud",
        y = "Latitud",
        title ="Cantidad de intervenciones con puestos móviles por comunas",  
        subtitle = "Ciudad Autónoma de Buenos Aires", 
        fill = "Número declarado de intervenciones",
        caption = "Fuente: IZLP, Mascotas de la Ciudad")

En el año 2019 hubo una fuerte presencia de los efectores públicos en la zona sur, con una distribución que recuerda a un Pareto, donde tres comunas, la 8(sur), seguida de la 4(sur) y 7(centro) concentran la mayoría de las intervenciones por parte de los efectores públicos con sus puestos móviles.

Se me ocurrió pensar que tal vez había alguna correlación entre la presencia de los servicios y los porcentajes de Necesidades Básicas Insatisfechas de las comunas, por eso los incluí en la base de comunas. Me parecía lógico que los servicios se concentren donde hay déficits en la CABA, pero a priori haciendo un gráfico de dispersión no hay tal relación. No ahonde lo suficiente como para encontrar algún índice que sean aptos para correlacionar más allá de lo planteado en el cálculo de REV.
Tampoco sé si un alto porcentaje de NBI por si solo sería un buen indicador de las comunas a las que habrían que darle prioridad de atención.

Mapas dinámicos:

Dada la magnitud de datos a representar y las múltiples aristas que se le pueden dar a medida que vamos analizando características de los servicios me parece una forma didáctica de mostrar la información en un mapa dinámico donde se puede hacer clic y obtener detalles, agregarle capas, etc. Estoy aprendiendo todavía así que están en construcción.

ciudad <- "caba" #lo que ponemos en el buscador

posicion <- 2 # La posición en que aparece entre los resultados de Nominatim

url <- URLencode(paste0("https://nominatim.openstreetmap.org/search.php?q=",
                        ciudad, "&polygon_geojson=1&format=geojson")) # sitio ciudad y archivo

destfile = paste0(tempdir(), "/limits.json") #archivo temporal donde se almacena lo que vamos a bajar

download.file(url, destfile) #pedido de descarga

limites <- st_read(destfile, quiet = TRUE)[posicion,] # los limites del mapa

MAPA DE TODOS LOS SERVICIOS: click para mas información

#Mapa con todos los puntos que prestan atencion clinica y resaltados en rojos los 24hs, queda super cargado cuando lo vez por primera vez peo el chiste de leaflet es que podes hacer zoom, buscar la zona donde estas y ver que opciones tenes disponibles:
mapa <- leaflet(limites) %>%
    addProviderTiles(providers$CartoDB.Positron) %>% #diseño blanco
    addPolygons(fill = NA) %>% 
    addMarkers(data = servet, popup = ~paste( "<h3>Efector:", unidad, "<br></h3>",
                              "Dirección:", direc, "<br>",
                              "<h3>Servicios:", servicio, "<br></h3>",
                              "<h3>Frecuencia:", periodicidad, "<br></h3>",
                              sep = " " )) %>% 
  addCircleMarkers(data = privadas_24hs,
                   color = "red",
              popup = ~paste( "<h3>Efector:", unidad, "<br></h3>",
                              "Dirección:", direc, "<br>",
                              "<h3>Servicios:", servicio, "<br></h3>",
                              "<h3>Frecuencia:", periodicidad, "<br></h3>",
                              sep = " " ))
mapa

MAPA DE LAS 24HS DISPONIBLES EN CABA: click para mas información

#Algunos tipos de mapas que me gustaron y se pueden usar cambiando el parametro:
#providers$OpenStreetMap.HOT #buen diseño
#Esri.WorldStreetMap #OldSchool

#Mapa de solamente los 24hs:
v24 <- leaflet(limites) %>%
addProviderTiles(providers$HikeBike.HikeBike) %>% #Detallado
  addPolygons(fill = NA) %>% 
  addCircleMarkers(data = privadas_24hs,
                   color = "blue",
              popup = ~paste( "<h3>Efector:", unidad, "<br></h3>",
                              "Dirección:", direc, "<br>",
                              "<h3>Servicios:", servicio, "<br></h3>",
                              "<h3>Frecuencia:", periodicidad, "<br></h3>",
                              sep = " " ))

v24

Para relacionar los datos que relevamos con la población que aspiran a satifacer, encontré una forma gráfica muy interesante para establecer la porción de territorio servida por cada punto de servicio, se puede dividir la superficie de la ciudad en “celdas de Voronoi”. La llamada “partición de Voronoi” es la subdivisión de un plano de acuerdo a la ubicación de un conjunto de sitios, donde a cada punto le corresponde una y sólo una celda. Uno de los atributos de la partición de Voronoi es que cada celda abarca los sitios del plano que se encuentran más cercanos a un punto particular y no a cualquier otro.
La metodología es utilizada con frecuencia para establecer áreas de servicio o cobertura y si pudiera afinarlo, lo cruzo con la info de la EAH para que estime la población de animales (con sus índices de no castrados, no vacunados, no desparasitadaos) y nos diría el número teórico de población animal que debería cubrir cada puesto.

En resumen, al calcular las celdas de Voronoi en torno a dos servicios, se estable el área donde residen la potencial poblacion objetivo, asumiendo que los mismos recurrirán al punto de venta más cercano.

MAPA DE LOS PUNTOS FRECUENTES CON AREA ESTIMADA DE COBERTURA: click para mas información

#Para finalizar algo que me quedo a mitad de camino pero es una idea que me parece genial para tener un parámetro que me diga que puntos de los servicios públicos hay que reforzar es el cálculo de las celdas de voronoi que hizo Antonio:

vetes_presencia <- servet %>%
  filter(!is.na(lon), !is.na(lat)) %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326)  %>% 
  # recortamos los puntos para retener sólo los que caen dentro de los límites de la ciudad
  filter(st_within(., limites, sparse = FALSE))

vetes_publicas <- puestos_frecuentes_geo %>% #nos basamos en los puestos PUBLICOS FRECUENTES
  mutate(Nro = row_number()) #esto lo hizo para hacer un join en el trabajo original, cosa que todavia no logré hacer en el mio 

celdas_voronoi <- vetes_publicas %>%
  st_union() %>% 
  st_voronoi()  %>% 
  st_cast() %>% 
  st_sf() %>% 
  # agegamos el identificador
  st_join(vetes_publicas["Nro"]) %>% #este es el espacial join que todavia no puedo hacer con mis datos
   st_intersection(st_geometry(limites)) # recortamos el resultado con la forma de los límites de la ciudad

celdas_voronoi

## Simple feature collection with 33 features and 1 field
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -58.53145 ymin: -34.70564 xmax: -58.33512 ymax: -34.52655
## geographic CRS: WGS 84
## First 10 features:
##    Nro                       geometry
## 1   18 POLYGON ((-58.4872 -34.5434...
## 2   15 POLYGON ((-58.52498 -34.601...
## 3   33 POLYGON ((-58.48891 -34.684...
## 4   20 POLYGON ((-58.44531 -34.535...
## 5   23 POLYGON ((-58.46115 -34.594...
## 6   31 POLYGON ((-58.46242 -34.589...
## 7   19 POLYGON ((-58.51239 -34.574...
## 8   21 POLYGON ((-58.47201 -34.576...
## 9   17 POLYGON ((-58.46003 -34.605...
## 10  11 POLYGON ((-58.50357 -34.621...

leaflet(limites) %>% 
  addProviderTiles(providers$HikeBike.HikeBike) %>%
  addPolygons(data = celdas_voronoi, fill = FALSE, color = "green") %>% 
  addMarkers(data = vetes_publicas, 
             popup = ~paste( "<h3>Efector:", unidad, "<br></h3>",
                              "Dirección:", direc, "<br>",
                              "<h3>Servicios:", servicio, "<br></h3>",
                              "<h3>Frecuencia:", periodicidad, "<br></h3>",
                              sep = " " )) 

Acá me quede, el siguiente paso en mi esquema mental seria calcular la población animal que cada punto tiene que atender para ver qué puesto de atención requiere refuerzo. Sin tener el dato tan exacto de población animal, pero sabiendo como está distribuida por el mapa que confeccionamos, es evidente que hay zonas donde un puesto único tiene una extensa área territorial con una considerable población animal que satisfacer… ¡Y recordemos que de estos puntos solo 4 tienen atención permanente! El resto es semanal, quincenal o mensual en puestos móviles que tienen muchas carencias propias de su poca infraestructura y horario reducido, sabiendo que los recursos que tenemos son acotados, una mejor gestión y organización de los elementos utilizables basado en la mejor evidencia disponible van a tender a mejorar los servicios ofrecidos a la población.
Estuve pesadeando a los que hicieron la EAH y me pasaron la base ampliada donde se pueden calcular los porcentajes por comuna, que es mas precio que Zona como lo publicaron ellos. Queda pendiente sacar 3 porcentajes particulares que podemos calcular sobre los datos que tomaron en las encuestas:
1. Animales no castrados.
2. Animales que no recibieron vacuna antirrábica en el último año.
3. Animales que no recibieron desparasitación en el ultimo año.
La idea es calcular un promedio de toda CABA (o de las zonas que estipulo la encuesta, no sé, todo es charlable) y ver cuáles comunas se desvian de esa media.
El detalle es que se pueden tomar las medidas como poblacion animal total o aprovechar el detalle que la EAH diferencia en caninos/felinos para lo cúal sabriamos más precisamente en que grupo tenemos que concentrarnos. Incluso está el dato de grupo etario al que pertenece, pero no me parecio relevante como para considerarlo en este analisis en particular.
Falta mas laburo para sacarle mayor provecho a los datos, pero en líneas generales creo que va encaminado, está todo abierto a sugerencias, gracias por leer todo jaja. ESPERO SU FEEDBACK!!